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We discuss the dynamics of a quantum phase transition in a spin-1 Bose-Einstein condensate when 
it is driven from the magnetized broken-symmetry phase to the unmagnetized "symmetric" polar 
phase. We determine where the condensate goes out of equilibrium as it approaches the critical point, 
and compute the condensate magnetization at the critical point. This is done within a quantum 
OO , Kibble-Zurek scheme traditionally employed in the context of symmetry-breaking quantum phase 

c7^ ■ transitions. Then we study the influence of the nonequilibrium dynamics near a critical point on 

^^ ' the condensate magnetization. In particular, when the quench stops at the critical point, nonlinear 

Cn \ oscillations of magnetization occur. They are characterized by a period and an amplitude that 

^.^. are inversely proportional. If we keep driving the condensate far away from the critical point 

j^ ' through the unmagnetized "symmetric" polar phase, the amplitude of magnetization oscillations 

slowly decreases reaching a non-zero asymptotic value. That process is described by the equation 
that can be mapped onto the classical mechanical problem of a particle moving under the influence 
of harmonic and "anti-friction" forces whose interplay leads to surprisingly simple fixed-amplitude 
oscillations. We obtain several scaling results relating the condensate magnetization to the quench 
rate, and verify numerically all analytical predictions. 
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I. INTRODUCTION 



C^ . Any symmetry breaking phase transition can be traversed in two opposite directions. The usual way is to start 

in the symmetric phase, and move into the phase with a broken symmetry. When this process happens on a finite 
timescale, critical slowing down will generally lead to the random local choices of the broken-symmetry vacua. This 

'^ ' in turn can result in the formation of topological defects [l|, [J- They will appear with the density set by the size 
5 ! of the regions that choose the same broken-symmetry vacuum. That size - and, hence, the density of the resulting 
non-equilibrium structures - can be deduced from the critical scalings of the relaxation time and the healing length. 
The obvious question that has not been addressed to date concerns traversing the transition from the phase in 
which the symmetry is broken to where it is restored. Such transition can be also induced on a finite timescale, which 

^^ ' will again bring into the discussion the critical scalings. Now, however, they will determine remaining excitation of 
the system rather than the size of broken-symmetry regions. 

We study here the dynamics of quantum phase transitions [3|, |j, |5|, la| rather than their classical (thermodynamical) 
counterparts that have received extensive theoretical [l|, Q and experimental 0, @] attention (see @ for a recent 

fT^ ' review). More precisely, we are interested in the dynamics of a ferromagnetic Bose-Einstein condensate (BEC) 
composed of spin-1 atoms [10|. Such a condensate is studied experimentally by several groups [UlllJ]- Similar issues 
may also arise in an even more complex case of higher spin condensates [l3|. The research on symmetry breaking 

J — I in the spin-1 BEC has been the subject of the seminal paper of the Berkeley group [11|. There the condensate 

C^ , was rapidly driven from the polar (symmetric) to the ferromagnetic (broken symmetry) phase by the decrease of 
a magnetic field. Formation of top olog ical defects (e.g., vortices) was observed. This work has triggered several 
theoretical investigations |ljj [la M, [l3, Ej [3 ■ ^^i particular, vortex density after a slow quench was found through 

^ ■ the extensive numerical simulations by Saito, Kawaguchi and Ueda [l3| to follow the scaling result obtained from the 
Jh \ quantum Kibble-Zurek mechanism [17] and the studies of the magnetization correlation functions [15l.ll§|. 
. - - I Below we consider dynamics of a ferromagnetic condensate driven through a transition "in the opposite direction" : 

from the broken-symmetry phase to the polar phase. As in [17| . we are interested in slow transitions that exhibit 
adiabatic-impulse behavior so that the condensate goes out of equilibrium close to the critical point. In this case the 
condensate excitation in the polar phase reflects the scalings of the critical regime. This scenario is supported by the 
following general discussion. Suppose that we change some dimensionless parameter q to drive the system towards 
the critical point at qc- Close to the phase boundary, the gap A in the excitation spectrum behaves as Ao((?c — qY" , 
where z and v are critical exponents [20|, q < Qc, and Ag is a constant with dimension of energy. As long as the 
system is driven far enough from the critical point, the gap is large and the system excitation does not occur. The 
relevant "reaction time" for the system is given by 

1-^ (1) 

The system undergoes adiabatic evolution when r is small compared to the timescale on which changes occur in the 
Hamiltonian. That timescale, in turn, is characterized by how fast the "instantaneous" excitation gap changes, or in 



other words by how fast the system is driven. It is simply given by 

A 



I rfA l 
I dt I 



(2) 



having a proper dimension of time. Comparison of the timescales ([T]) and ([2|), 

h A 



A i^r 

I dt I 

gives the location of the border between adiabatic - non-adiabatic behavior at some q ^ q^ — q. Assuming, that 
■^q{t) = Tq^ near a critical point {tq^ is the speed of transition between the phases) we get that 

The qualitative picture of system dynamics is as follows. For qc — q{t) > q the system undergoes adiabatic evolution 
following instantaneous changes to its Hamiltonian. When qc — q(t) ^ q the evolution ceases to be adiabatic and to 
a first approximation one can expect the impulse stage, where the state of the system does not change - remains 
"frozen". As a result, system properties near the critical point are determined by its ground state at q — q^, ^ q, 
where q is computed from the quench time tq and the product of the critical exponents: Eq. ([3]). This approach 
was already successfully used to study the Landau-Zener dynamics [5|, the dynamics of the quantum Ising model 
[31, and the ferromagnetic condensate dynamics during the symmetry-breaking transitions from the polar to the 
broken-symmetry phase [17| . 

We begin our study in the next section, where we discuss the basics of a quantum phase transition in the mean-field 
model that we consider. Section IIIII presents how the condensate gets excited when driven through the broken- 
symmetry phase. Section IIVI is devoted to studies of magnetization oscillations after quench. We point out the 
observables that contain the information about condensate excitation at the phase boundary. Section |V] determines 
the range of applicability of the Single Mode Approximation (SMA) used in Sections IIIII and IIVI Finally, Section IVII 
provides the summary of the paper. 

II. THE MODEL 

The energy of the ferromagnetic condensate placed in an external, homogeneous, magnetic field aligned in the z 
direction is given in the mean-field approximation as 

£m=fdr^\V^\^ + ^{^^^y + Q{^\F!\^) + 'j- Yl (*l^"l*>'> (4) 

a—x^y,z 

where Fa is the standard spin-1 matrix, and co > and ci < provide the strength of spin independent and 
spin-dependent interactions, respectively: 
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with as being the s-wave scattering length in the total spin S channel, M is the atom mass, while Q (proportional 
to square of magnetic field imposed on the condensate) is the prefactor of the quadratic Zeeman shift resulting from 
atom-magnetic field interactions. The linear Zeeman term is skipped in Eq. ^ because it has a trivial effect on 
condensate dynamics. Indeed, it only rotates the condensate magnetization around the z axis with the Larmor 
frequency. This rotation, of course, can be removed from theoretical studies by going to a reference frame rotating 
with the Larmor frequency (Sec. II. A of Ref. p^). 

The wave function has three condensate components, ipm, corresponding to ?ti = 0,±1 projections of spin-1 onto 
the magnetic field: 



* = V^o , / dr*^* = N, (5) 




where N is the total number of atoms. In subsequent calculations we will use the phases Xm defined as: ipm = 

|V'm|exp(«Xm)- 

The simplest meaningful approach to the dynamics of a ferromagnetic Bose-Einstein condensate is provided by the 
Single Mode Approximation where the translationally invariant mean-field approach is used. This simplification, while 
obviously inapplicable for symmetry-breaking phase transitions where the orientation of magnetization is explicitly 
position-dependent [ITI, Ii3| , can be successfully used for the transitions considered here where the system starts from 
a uniform broken-symmetry ground state (GS). The influence of inhomogeneities will be illustrated in Sec. FVl where 
we will compare simulations done with the SMA to full mean- field calculations. In the SMA both the gradient term 
and the density-density interaction term ~ cq are irrelevant and so are skipped in Sees. [TIllIVI These terms, however, 
reappear in the numerical calculations in Sec. El 

As the ferromagnetic condensate is described in the SMA by translationally invariant spinor ([5]), it is instructive to 
parametrize the coupling constant Q from ^ by 

Q = qn\ci\, 

where n = ^^'I' is the atom density, and g is a dimensionless parameter proportional to the square of the magnetic 
field. This parameter will be changed to drive the condensate from one quantum phase to another. 
The system longitudinal (fz) and transverse (f^ and fy) magnetizations read 

/a = (*|i^a|*), a = x,y,z. 
We consider the fz — case, i.e., |-0i| = |'0-i|- This condition is dynamically conserved during all the evolutions 



considered in Sees. Hill and HVl because -^ J drfz{r,t) = 0, which in the SMA reduces to ^fz{t) = 0. The /; 
components are conveniently combined to a complex transverse magnetization 



lx,y 



fT = !x+ ify = V2(V'iV'o + V'-iV'o) = 2%/2|?/'i||Vo| cos ( xo 



Xi+X-1 \ t4, 



X-1 -Xi 



(6) 



The quantum phases of the ferromagnetic condensate were recently discussed in [2l| . The condensate is in the polar 
phase when q > 2. There the ground state is given by 




so that fx = fy = 0. For < g < 2 the system is in the broken-symmetry phase where the GS spinor reads 
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and 



1 - ^e't 



It = n^ 
The system dynamics in the SMA is governed by [20] 

thi^i = ci|Vo|Vi+ciVoVi*e-''^4 

</) = 0, 



(8) 



n\ci\q-)pi 



(9) 



where d/dt is denoted with a "dot" and we have used the identity 1-011 = IV"-!! to arrive at Eqs. ^. The initial 
condition for time evolutions in Sees. IIIII and lTyi is provided by the spinor state (O taken at q = 0. The last equation 
of ^ states that the orientation of the transverse magnetization on the {x, y) plane is conserved throughout the 
evolution and so it is fixed by the initial conditions. This constraint prohibits full restoration of the symmetry in the 
state evolved from the broken-symmetry phase to the symmetric polar phase. In particular. 



/:r(t)sin(</)) = fy{t)coii{(j)). 



We will discuss below the dynamics driven by the variation of q that can be achieved with a proper manipulation of 
the magnetic field imposed on the condensate. A linear increase of the magnetic field strength in time is well within 
the experimental capabilities [ll[. Due to the quadratic Zeeman coupling, it results in 

"'>4(4)^^ (10) 

This is the time dependence we will assume in all subsequent calculations. Close to the phase transition point, where 
t = Atq — 5t, we have 

For the transitions considered here the system will cease to adiabatically follow the ramp up of q near a critical point 
where the above linearization holds: the familiar time-dependence from the Kibble-Zurek theory emerges near the 
phase boundary with Tq providing the quench rate j23|], or in other words, the speed in the parameter space of 
driving the system through the critical point. 

III. SYSTEM EXCITATION 

For small perturbations around the GS of the broken-symmetry phase one finds three Bogolubov modes as in [2l| : 
two gapless modes and one gapped. We do not consider analytically the dynamics induced by the excitation of 
the gapless modes. Instead, we focus on the gapped mode and show that its excitation is responsible for changes 
of the magnitude of the system transverse magnetization during non-eq uilibrium dynamics. This mode in the long 
wavelength limit, obviously relevant for the SMA, has the energy gap [2l| 

A-2n|ci|Vl-gV4- (11) 

As we drive the system to the critical point changing q{t) we expect qualitatively the following dynamics. First, for 
q{t) small enough the evolution is adiabatic with respect to the gapped mode - far enough from the critical point the 
energy cost of exciting the gapped mode is too large. However, as the critical point is approached, the gap becomes 
too small to allow for adiabatic evolution and the system starts to populate the gapped mode: the non-adiabatic 
dynamics starts. Below we assume that the separation of the two regimes takes place at q = 2 — q, and determine the 
location of q as a function of the quench rate. This can be done as in Sec. HI Comparison of the timescales ([T]) and 
^ brings us to the (approximate) equation 



A 



A(g) |A| ^ 
This equation is illustrated in Fig. [1] Its solution in the slow transition limit, g ^ 1, gives 



(12) 



1 / \ 2/3 
1 / To ^ ' 



42/3 VtqJ ' ^"- n\c^\- ^^^^ 

This, of course, corresponds to ^ with zi' = 1/2 and Aq — 2n\ci\. A simple estimation of the constant tq by the 
peak density in the Berkeley experiment, 2.8 x 10^'^cm~^, gives tq w 16ms. In this section and Sec. IIVI all the results 
can be obtained without assuming any particular value of tq . 

Qualitatively, prediction (|13p is in agreement with what we expect: the slower the quench is (the larger tq is), the 
closer to the critical point the condensate can be driven adiabatically. 

Quantitatively, solution p3p turns out to be remarkably accurate in the wide range of quench times tq as depicted 
in Fig. [2l In that plot we define the instant q to be the distance from a critical point when the departure of system 
transverse magnetization from a static (ground state) prediction starts exceeding some threshold, e.g., 5% of the 
instantaneous GS value. From the fit there we get that 
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FIG. 1: (color) Illustration of Eq. p2|) in units of Tq (|13p on the vertical axis. The plot is for tq/tq = 1. 
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FIG. 2: (color) Departures of the transverse magnetization from the static prediction during evolution in the broken-symmetry 
phase, l/rlstatic refers to the squared modulus of ((S)). We arbitrarily define that when these departures exceed 5%, the system 
leaves the adiabatic regime and starts the non-equilibrium dynamics (other reasonable thresholds < 10% give the same scaling 
result). The inset shows the result of power law fit pointing to g = 1.05 x (to/tq)"'^^^*'''^''^ , which compares well to theoretical 
prediction (|13p . The fit is presented as a solid line, while numerics comes as pluses. Note the log- log scale on the inset. The fit 
was done for tq/tq = 30 • • • 200. The larger tq's are taken to the fit, the closer scaling exponent approaches 2/3. On the main 
plot we have a simulation for tq/tq = 30. 



in excellent agreement with (J13p with respect to the scaling exponent. The prefactor, which obviously depends on 
the threshold used for the numerical determination of q is of the same order as in (I13p . Additional details are seen in 
Fig. [21 The error in the determination of the scaling exponent is one standard deviation coming from a linear fit on 
a log-log plot. All fitting errors in this paper are determined in this way. 

Let's look at the system magnetization at the critical point. To proceed further it is convenient to define the 
condensate energy density as 



E- 



-^I/tP 



2n\ci\q\ipi 



(14) 



which is equal in the SMA to the contribution of the last two terms in Eq. ^ to the condensate energy per volume 
(the first two terms in Eq. ([4]) are unimportant in the SMA as the gradient term equals zero while the density-density 
interaction term ^ co is a constant). In Fig. [3] we observe that |/tP ~ Tq and E ^ Tq . Putting these two 
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FIG. 3: (color) Scaling of the condensate properties at the critical point. The upper data presents |/t| /n , where /t is the 
transverse magnetization, while the lower data is for E/n^\ci\, where E (energy density) is given by (|14|) . The pluses come 
from numerics, while the solid lines present the curves ~ (to/tq)'' with 7 = 2/3 (4/3) for the upper (lower) data. This confirms 
that at the critical point 1/tP 
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scalings into Eq. p^ we see that 
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at the critical point in the leading, i.e., Tq , order in the quench rate. 

It is now interesting to ask whether there is any relation between the condensate magnetization at the point where 
the system goes out of equilibrium, q ~ 2 — q, and the condensate magnetization at the critical point, q = 2. The 
condensate magnetization at q = 2 — q equals (in the slow transition limit where ^ <C 1) 



\fT\Vn' = l-{2-q)y4: 



-2/3 



so Ifriq — 2 — <7)P scales with the quench rate in the same way as Ifriq — 2) p. This resembles the adiabatic-impulse 
simplification of dynamics [3|, |5| , where it is assumed that the evolution of a system undergoing a phase transition 
is either adiabatic or impulse, and the impulse stage implies no changes in the state of the system. Here this simple 
picture gives a correct scaling of the system magnetization at the critical point from the value of q - the location 
where the system leaves the adiabatic regime. This is a good approximation despite the fact that the condensate does 
not enter an ideal "impulse" regime. Indeed, the wave-function is not "frozen" and does change from q — 2 — q to 
(7 = 2. Remarkably, however, these changes are such that instead of getting magnetization equal to |/t('Z = 2 — g)p 
at the critical point (ideal impulse dynamics) we get const ■ |/t(9 = 2 — q)P, where the constant is TQ-independent. 
This character of condensate evolution toward a critical point is universal with respect to quench time. It is depicted 
in Fig. S 

Now we would like to address how these critical scalings can be experimentally observed. One obvious way is to 
watch the condensate dynamics at different times as was done in the Berkeley experiment |ll| . Another approach 
involving interesting physics would be to stop a change of q and watch subsequent magnetization dynamics. This 
approach will be discussed in the next section. 



IV. DYNAMICS OF MAGNETIZATION AFTER THE QUENCH 



Suppose that one stops driving the system at some q > 2 and lets it evolve freely, which leads to magnetization 
oscillations. We would like to investigate here whether the information about the system excitation at the critical 
point can be extracted from the amplitude and/or period of these oscillations. In other words, we want to find out 
what are the signatures of the nonequilibrium dynamics in the broken-symmetry phase that may be observable in the 
polar phase. This problem is of both experimental and theoretical interest. 
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FIG. 4: (color) Illustration of the scaling of departures of the condensate transverse magnetization, |/tP, from the static 
prediction, l/rlatatio given by the squared modulus of (|8]). The solid black line is for tq/tq = 30, while red dots show data for 
■'"q/'^o = 60. The plot shows that the departures of the transverse magnetization from a GS prediction near a critical point 
scale as Tq after rescaling of the distance from the critical point, 2 — g, by a Tq factor. 




FIG. 5: Time evolution of the magnetization during the quench that stops at the critical point. This plot is for tq/tq = 10. 
The critical point is reached at t/ro — Atq/tq — 40, so q{t < 40ro) is given by (|10|l while q{t > 40ro) — 2. The system excitation 
due to approaching the critical point is visible from the nonzero magnetization at t/ro > 40: if the evolution would be adiabatic, 
the condensate magnetization at the critical point would be zero and no magnetization oscillations would be present on the 
plot. 



A. Quench stops at the critical point 



Suppose the condensate is driven from q ~ to the critical point with (jlOp . The critical point is reached at t — Atq 
and then q{t > Atq) — 2, i.e., free (without any driving) evolution takes place. The evolution of magnetization for 
such a problem is presented in Fig. [51 where periodic oscillations are easily seen. Our aim is to describe them and 
show that both their amplitude and their period contain information on how the condensate was excited at the critical 
point. 

In this case one can derive the following equation for the evolution of the magnetization: 



frit) 
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E + 2 
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ft2 



■I/tWH/tW, 



(15) 



where E is the system energy density (J14p at the critical point (in fact, also at any other time since lack of parameter 



8 

changes results in energy conservation). This equation can be obtained straightforwardly from Eqs. ([6]) and ([9]). It 
is an interesting equation: in the limit of slow transition, when the amplitude of oscillations becomes very small, the 
nonlinear term dominates the physics rather than being negligible compared to the ^ E term promoting harmonic 
oscillations. Indeed, from Sec. Illll we see that |ci|_E/|ciP|/tP ~ {'''o/tq)^^^ ^ 1 for slow enough quench. 
The exact solution reads 



/T(t)=Acn(^i + 0,«; 



(16) 



where en is the Jacobi cosine 124 



1 



l + aV2' 



o^iMP 



and A and 9 are given by initial conditions. They both might be found from a fit to experimental data. In particular, 
|Ap is given by the amplitude of free |/tP oscillations, because the Jacobi cosine ([T6| oscillates periodically between 
±1. 

The Jacobi cosine was formerly used in the spin-1 condensate context in [12l |25|. It is periodic - cn(a;, k) = 
cn(a; + 4:K{k), k) ~ where K{k) is the complete elliptic integral of the first kind. In the time domain this periodicity 
translates into 



At = 4k 



Aci 



K{k) 



(17) 



To simplify this expression we note that from the fact that the Jacobi cosine in (1161) is bounded by ±1, Jt and A are 



of the same order and so we expect \A\ ^ 
transitions) a ^ (tq/tq)^^^ '3> 1 and so k 



^-1/3 



Moreover, as E 



-4/3 



we get that in the limit of large tq (slow 



l/\/2- This last result allows us to write 



At 



r(i/4)= 



Aci 



(18) 



This finding is quite interesting: the periodicity is amplitude-dependent unlike in the harmonic oscillator case. This 
simple result predicting the inverse proportionality of the period and the amplitude of oscillations is accurate to about 
1% for experimentally relevant quenches as depicted in inset (a) of Fig. [6l 

The numerical results on the oscillation period At when the quench stops at g = 2 are presented in Fig. [S) From 

the fit we see that At scales as tJ in accordance with expression (fT5|) supplemented with the above observation that 



\A\ 



-1/3 



(proven numerically in the inset (b) of Fig. [5]) . 



B. Quench stops after passing the critical point 



According to our numerical simulations there are three stages of such evolution depicted in Fig. [Tj (i) the system is 
driven to the critical point through the entire broken symmetry phase and its magnetization at the critical point scales 
as I/tP ~ Tq (only barely noticeable magnetization oscillations are present - see Ref. [23[); (ii) the condensate is 
driven by a change of q, Eq. (jlOp . in the polar phase and we observe there damped magnetization oscillations; (iii) 
the free evolution after we stop ramping up q{t) takes place and periodic magnetization oscillations appear. These 
oscillations are described by 



/t = - 



r,'q{q-2)+A^-^E- 
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(19) 



where again E ^Q and Eqs. ([6]) and ([9]) have been employed in derivation of ([19 
Jacobi cosine function exists and has form (J16p with 



The exact solution in terms of the 



q{q~2)h^/2T^ + 2E\ci\ 



Now, in the limit of fixed q > 2 and large tq (slow transition) we have a different behavior of the Jacobi cosine than 
above. Namely, now a ^ (tq/tq)^^^ <C 1 which results in k ^ and cn(»,K) -^ cos(»). In this limit the last two 



38 



28- 



< 



18 







' ' 


' ' 


1 


/ 


1% 


\ 


1 ' 1 ' 1 


1 ' 1 ' 

(a)_ 


0.25% 


- 


I.I.I 


,1 






30 


60 90 120 


150 180^ 




^q/^0 


^ 


y" 






0.25 


+ 


, 1 I 1 1 1 


(b)'- 


^ 


lAI/n 


0.15 




1 1 . 1 1 1 


1 


X 






30 


60 100 

1 'q^'o 


200 



30 



60 , 100 

T /T 
Q 



200 



FIG. 6: (color) Here we show the properties of the system when the change of q stops at the critical point and then there is 
a free evolution going on. Note that there is a logarithmic scale on the vertical and horizontal axes of the main plot and the 
inset (b). Main plot: the periodicity interval Ai of fx,y{t) during free oscillations. The fit (red line) to numerics (pluses) gives 
At ~ .j-0.3350io.oooi jjjgg^ (a,): the difference between (|T5|| and (fT7)) divided by ([T7)). i.e., relative departures of the approximate 
formula for oscillation period from the exact result. The approximate formula (|18p overestimates the exact result by about one 

ich time tq . The fit (green line) to 
in agreement with adiabatic - impulse simplification of the condensate dynamics. 



percent in a wide range of tq/tq's. Inset (b): dependence of \A\ from (|16p on quench time tq. The fit (green line) to numerics 
(pluses) gives \A\ ~ t-q°-3304^o.oooi 



terms in (J19p can be neglected and the Jacobi cosine turns into a normal cosine function: harmonic oscillations show 
up with a repetition period 



At 



2ttto 



^q{q - 2) 



(20) 



Naturally, the transition between Jacobi cosine oscillations and typical harmonic dynamics is gradual and can be 
traced quantitatively with the exact solution p6)) . We can compare ([20]) to the numerics. For the evolutions where 
the increase of q stops at g = 3, as is the case in Fig. [71 we get from ([20]) that At = 2 x 1.814to, while from numerics 
we obtain 2 x 1.804to for tq/tq = 30 and 2 x 1.812ro for tq/tq = 200 (notice that |/tP oscillates with 1/2 of the 
period ([20)1 and so our numerical results extracted from |/tP oscillations were multiplied by a factor of two). As 
we see, the larger tq, the better the agreement because the nonlinear corrections become smaller. This is in fact 
bad news: as pop is rg-independent, we can no longer use repetition period At to investigate the dynamics of the 
quantum phase transition. 

Fortunately, however, behavior of the amplitude of free magnetizations oscillations provides an easily visible signa- 
ture of the non-equilibrium dynamics. As shown in the inset of Fig. [3 the amplitude of |/tP scales as Tq . This 
scaling results from two observations. First, the driving in the polar phase starting from q ^ 2 and proceeding to 
g > 3 damps the amplitude of |/tP by a factor of Tq . We have verified this numerically by simulating different 
evolutions that begin from a fixed, tq independent, state. Second, the scaling of |/tP at the phase boundary is Tq . 

Combining these two results one easily justifies Tq scaling from the numerics. 

Finally, we comment on what happens when the system is continuously driven through the polar phase. Combining 
Eqs. ([nj and ^ one easily arrives at 



h 



''q{t){q{t)-2)h 



^ q{t) 
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(21) 



that has to be solved with initial conditions given at instant ti by friti) and friti) = iqV2Tf^^ (tpl^po — ^-i'^o)\ti- 
Two remarks are in order. First, the same initial conditions apply to (jlSp and (|19p . Second, equation (|2ip is no longer 
self-consistent: an additional equation for ijji{t) dynamics has to be solved simultaneously unless the term ~ 1-01(^)1^ 
is negligible. 

As we have observed in Fig. [7] the magnetization oscillations arc damped when the system is driven in the polar 
phase. It turns out, however, that when we continue driving the condensate, the amplitude of magnetization oscilla- 
tions reaches some nonzero asymptotic value, scaling as Tq , instead of decreasing to zero. This is supported by Fig. 
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FIG. 7: (color) Dynamics of the system driven by the change of q{t) given by (|10p until the point q — 3, and then undergoing 
a free evolution with q time independent. The main plot shows data for tq/tq — 30. The free evolution starts at t/ro = 
\/24rQ/ro « 147. Different stages of evolution are encoded in colors: black (evolution in broken-symmetry phase), red 
(dynamics in polar phase induced by change of q), green (free dynamics in the polar phase, q = 0). Inset: the scaling of the 
amplitude of [frit)]'^ oscillations during free evolution at g = 3 (depicted in green). Pluses come from numerics, while the line 
is the fit corresponding to ~ .^-looisio.oooi gcaijng 



0.025 




FIG. 8: Maxima of the system magnetization, |/t| , during an uninterrupted driving through the polar phase with q{t) given 
by (|10p . The solid line comes from numerics and provides the upper envelope for oscillations of |/T(i)|^. The data is for one 
run with tq/tq — 30. The dashed line is the asymptotic (large q) value for the solid line. It is obtained from a linear fit 
(Maxima of 1./tP)/w'^ — a + b- for q ~ 100 data. The fit gives a « 0.0074 plotted as a dashed line. Additionally, we have 

verified numerically, that for different tq's the following scaling holds: a ~ Tq^ j which can be explained in the same way as 
the fitting result from the inset of Fig. [T] 



|8]and the following analysis. First, we assume that q{t) ^ 2 such that it is safe to approximate q{t){q{t) — 2) by q{t)'^. 
Second, we neglect the second term in (|2ip as it is small compared to the first one in the slow transition limit. After 
that we end up with the following equation that can be solved exactly 



/.^-vW/. + f/- 



(22) 



This is precisely a driven harmonic oscillator with anti-friction since q{t)/q{t) is positive in our case. In a classical 
mechanical picture we may imagine that (j22p describes a particle in a harmonic potential under the influence of a 
force that acts along its velocity. Such a force fights the tendency to squeeze particle motion due to increase of the 
frequency of the harmonic oscillator. Interestingly, it turns out that the particle perfectly maintains the amplitude of 
its motion because the exact solution of ([H]), constrained hy (j) = from ^, is 

/t = Ccos ( — f dtq{t) + e] , C = 0. 
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Additionally, we stress the fact that this solution is valid for any smooth q(t) dependence and not just (jlOp . In 
particular, in the opposite limit when q{t) decreases in time the "anti-friction" term works as a special friction term 
that slows the particle down so that it oscillates with constant amplitude instead of spreading out. 

This phenomenon can be qualitatively explained as follows. As it happens at q(t) ^ 2 we can safely assume that 
the magnon modes of the system do not become excited any more during driving. This is so because they have energy 
gap that becomes large far away from the critical point. These modes are responsible for the scattering of atoms 
between m — and m — ±1 condensate components [21|. In the absence of that process the populations of to = 0, ±1 
sublevels should be constant, which implies that 

|/Tp = 8|V'inVopcos2L-^^^;^ 



will undergo fixed amplitude oscillations induced by rotation of the condensate phases. 

V. DYNAMICS IN A SLIGHTLY INHOMOGENEOUS SYSTEM 

In this section we would like to compare the results obtained from a Single Mode Approximation enforcing the 
translational symmetry onto wave-function of the system to the more experimentally relevant inhomogeneous ( "dis- 
ordered" ) problem. To this aim we assume that the atom cloud is placed in a box-like trap and impose density and 
phase fluctuations onto it. These may come from experimental imperfections and quantum fluctuations. 

The results presented here come from numerics done in a one-dimensional configuration and so we present an 
explicitly ID description below unless stated otherwise. We consider untrapped system: atoms in the box as in the 
experiment [2^ done with spinless bosons. We estimate the parameters for these ID simulations as in [l3| taking into 
account the experimental setup of [ll|. Below L « 400/i7TT- is the size of the box keeping N ^ 2 x 10^ atoms confined 
along the z direction. 

We introduce disorder by modifying the initial wave function at q{t = 0) = in the following way: 

/ ^M^)e^''' \ i-L 

where the mode amplitude is proportional to: 







\ 






9<'"'^M^^..( ^^~^' 



(™)n2' 



1 + <5 > ' e^vr' rr> exp - ^ ; ' . (23) 



=1 



There S gives the strength of the disorder, T]n' , x±i & [0, 27r) are random phases, r„ G [0, 1] are random amplitudes, 
and Zn G [0, L] are random positions of gaussian disturbances. All these parameters are generated with uniform 
probability density. We have chosen a to be equal to the spin healing length ^s and assumed that the average spacing 
between Gaussian perturbations is 2^^- The spin healing length in the experiment |11| is 2.4y^TO, so that k, the number 
of Gaussians in ([23|) . is about 80 (L/2^s). The meaning of S is that it provides a characteristic relative-to-background 
size of density variations. A typical disordered pattern in ([^5]) is depicted in Fig. ^ Now our goal is to answer what 
is the range of the 6 parameter that allows for observation of qualitatively the same physics as in the homogeneous 
system described in Sees. IIIII and ITVl 

To look at the dynamics in the broken-symmetry phase, we have done calculations for different S G (0,0.2). For 
each S we randomly generate many-times the pattern (|23p and evolve the system to the critical point at different 
quench times tq . The patterns are independently generated for every evolution, which shall resemble the experimental 
situation where the structure of "imperfections" fluctuates from run to run. Then for each pair (5, tq) we average the 
mean system magnetization over Nr = 30 runs. 



^ N L 



where ff^ {z,t) is the condensate transverse magnetization in the i-th run. 

We look at two quantities: the distance qirq) from the critical point where the system starts the non-adiabatic 
evolution due to approaching the critical point, and the scaling of the condensate magnetization at the critical point. 
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FIG. 9: Plot of \(j>rn{z)\'^ for a typical realization of disorder ([23}. The parameter 5 = 10%. 





A 


a 


S = 0.5% 


0.65 ±0.01 


0.650 ± 0.001 


5 = 1% 


0.67 ±0.01 


0.647 ±0.001 


S = 2% 


0.65 ±0.01 


0.60 ± 0.02 


(5 = 3% 


0.64 ± 0.02 





TABLE I: The scaling exponents A, g ~ ''"o^i ^.nd a, {|/t|^) at the critical point ~ t^ 

deviation) come from a linear fit on a log-log plot to data in the range of tq — Tims ■ 

{S = 0) points to A, cr « 2/3. The definition of q is the same as in Fig. [2]except (|/tP) is used instead of J/t 



The values and errors (one standard 
1.1s. The Single Mode Approximation 



Both are defined in the same way as in Sec. IIIII except that now we use (|/tP) instead of |/tP- Our results suggest 
that the scaling exponents in the disordered system can match the SMA predictions within a few percent accuracy 
for S < 2%. Quantitatively, as we see in Table [H the scaling exponent A given hy q ^ Tq^ stays reasonably close 
to the 2/3 value for the homogeneous system for S as large as 3%. For larger S, e.g., 6 — 5% the numerical data 
clearly departs from the power law. The scaling exponent ct, defined at the critical point as (|/tP) ^ '^q^ ^ deviates 
noticeably from the SMA prediction already for 8 = 2% (Table [l|. For larger S, e.g., S = 3% numerics does not follow 
the power law anymore. 

There are at least two reasons for discrepancies between disordered and homogeneous results. First and most 
importantly, the presence of disorder perturbs magnetization affecting determination of both exponents from Table 
|T1 Even when these departures from a GS are relatively small compared to the GS magnetization when the system 
starts time evolution at g = 0, they can be significant when compared to the condensate magnetization close to or at 
the critical point where scaling exponents are determined. Second, our prediction of the scaling exponents is based on 
translationally invariant theory, i.e., momentum fc = problem, while the introduction of disorder leads to population 
of fc 7^ modes. A more involved description with the gap A pTJ) being a function of k should be used when the 
population oi k j^ modes becomes significant. 

Now we would like to look at the system dynamics in the polar phase. In Sec. IIVI we have observed that the 
transverse magnetization undergoes periodic oscillations in a homogeneous problem after stop of driving. Now we 
see dephasing dynamics spoiling this picture after a time inversely proportional to the disorder magnitude and the 
quench time tq. The latter observation results from the fact that the slower we go at fixed S, the more significant the 
disorder is compared to the system magnetization in the polar phase. 

Fig. [TO] presents the evolution with tq = 74ms and i5 = 1%. The change of q stops at the critical point there. For 
these parameters the scaling exponents are very close to the Single Mode Approximation result (see Table |T| . As is 
illustrated in Fig. [TOl however, dephasing takes place after about 7 — 9 oscillations. To observe more magnetization 
oscillations closely following the SMA result, for twice that time, S has to be about one order of magnitude smaller. 
Therefore, the system dynamics is quite sensitive to inhomogeneities when the condensate is left at the critical point. 
This should result from the coupling between different k modes due to nonlinear terms in the evolution equation. 

Fig. [11] presents the evolution where driving stops at q — 3 (similarly as in Fig. [7]). The parameters and the initial 
state taken for this simulation are the same as in Fig. 1101 This time the undriven dynamics takes place in the regime 
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FIG. 10; (color) Dynamics of the transverse magnetization in the inhomogeneous system during a quench with tq = 74ms. 
The change of q{t) stops at the critical point: q{t > irg ~ 0.3s) = 2. The thick black line is for S — 1%. The thin green line 
presents 5 = result, i.e., the Single Mode Approximation outcome. The inset shows long time behavior of the black thick 
curve from the main plot. The initial state for that evolution is the same as in Fig. 1111 




FIG. 11: (color) Dynamics of the transverse magnetization in the inhomogeneous system during a quench with tq = 74ms. 
The change of q{t) stops at 3 in the polar phase, q{t > \/24tq ~ 0.36s) — 3, and 5 — 1%. The insets show enlarged data 
from the main plot corresponding to the regions depicted with boxes. The period of oscillations is so small that individual 
oscillations are not resolved in the main plot and the upper inset. The lower inset does not show the SMA result as it is 
practically indistinguishable from the inhomogeneous calculation. Arrows indicate which inset belongs to what box. The initial 
state for that evolution is the same as in Fig. 1101 



where nonlinear couplings are small. As a result, the dephasing occurs after hundreds of oscillations even though the 
disorder amplitude, 5 — 1%, is the same as in Fig. [TU] 

To explain dephasing when the evolution stops away from a critical point, we can linearize the coupled Gross- 
Pitaevskii equations. Below we do it in an explicitly 3D configuration to use the same notation as in Sees. Illll and HVl 
(the ID result is qualitatively the same). We assume that 



(5Vi(r,t) 
5V^_i(r,t) 



exp 



-i' 



where the chemical potential is ^ = nco, \54>m\ ^ 1, and J dY{5'4>i^ + Sijj'^) = to keep J dr^^^ — N + 0{6ijjf-,^). In 
leading order 

/T = \/2n(<5Vi+5V'-i), 
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and its dynamics at fixed q is given by 

It = -^^'fT + {<!- 1)^VVt - r,-\iq - 2)/t, 
where the last term is known from P^ . That can be solved by the substitution fx ~ cos(kr — uj(k)t), giving 



hu;{k) = ^el + ek2{q - l)n\c,\ + n^\c,\^q{q ~ 2), 

where e/c — h^k"^ /2M . When disorder is present, different k modes will be occupied and they will oscillate at different 
frequencies aj(fc), which will result in the dephasing observed in Fig. 1111 The expression for Lo{k) shall not be used at 
q = 2, where linearized theory fails as was shown within the SMA in Sec. IIVI 

VI. SUMMARY 

We have analyzed transitions from the broken-symmetry phase to the polar phase occurring on the finite timescale 
TQ. Our focus was on slow transitions, when the condensate goes out of equilibrium near a critical point. We have 
presented a simple theory predicting that the evolution will cease to be adiabatic before reaching the polar phase at the 

distance ''^ Tq from the critical point. This result was found to be in excellent agreement with numerics. Applying 
the basic assumptions of the adiabatic-impulse approach 3, |5(, which originates from the Kibble-Zurek theory of 
nonequilibrium dynamics of classical phase transitions [l|, y], we have explained why the transverse magnetization 

2/3 

of the condensate driven to the critical point scales as Tq as well. Subsequently we have analyzed how the 
latter result can be experimentally extracted from the observation of magnetization oscillations. Three cases were 
considered. First, we have studied what happens when the quench stops at the critical point and the system follows 
free evolution. It was shown that the periodic oscillations appear with the repetition period coupled to the amplitude 

1/3 

of transverse magnetization oscillations. The repetition period was found to scale as Tq , where the exponent was 
directly related to the scaling of the transverse magnetization at the critical point. That nonlinear dynamics was 
exactly described within the mean-field approach in terms of Jacobi elliptic functions. Second, we have studied what 
happens when the system is driven into the polar phase, and then undergoes a free evolution. In that case we have 
shown that its excitation at the critical point can be easily extracted from the scaling of the amplitude of the free 
magnetization oscillations given by Tq . Third, we have studied what happens when the condensate is driven without 
any interruptions not only in the broken-symmetry phase, but also in the polar phase. In this case the amplitude of 
driven magnetization oscillations reaches a non-zero asymptotic value: the condensate magnetization is described as 
an unusual harmonic oscillator with anti-friction term that perfectly cancels the squeezing induced by the increase of 
the harmonic oscillator frequency. 

All the above results were obtained within a translationally invariant Single Mode Approximation whose range of 
applicability was determined by introducing a controlled disorder into an initial wave function. We have found out 
that the amount of disorder has to be quite small to have the homogeneous system predictions applicable. 

Future extensions of this work will include studies of quantum phase transitions dynamics in a harmonically trapped 
ferromagnetic condensate and investigations of the role of quantum fluctuations on the condensate dynamics. 
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